/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using System;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal abstract class AbstractSvd
    {
        protected int mColumns;
        private bool mComputed;
        protected bool mComputeVectors;
        protected bool mConverged;
        protected int mRank;
        protected int mRows;

        protected Vector mS;
        protected Matrix mU;
        protected Matrix mV;

        protected AbstractSvd(bool computeVectors)
        {
            mComputeVectors = computeVectors;
        }

        public bool Converged()
        {
            Compute();
            return mConverged;
        }

        public double Norm2()
        {
            return mS[0];
        }

        public double ConditionNumber()
        {
            int tmp = System.Math.Min(mRows, mColumns) - 1;
            return mS[0]/mS[tmp];
        }

        public int Rank()
        {
            return mRank;
        }

        public Matrix U()
        {
            return mComputeVectors ? mU.Clone() : null;
        }

        public Matrix VT()
        {
            return mComputeVectors ? mV.Clone() : null;
        }

        public Matrix W()
        {
            int rows = mU.Rows;
            int columns = mV.Columns;
            Matrix result = mU.CreateMatrix(rows, columns);
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    if (i == j)
                    {
                        result.ValueAt(i, i, mS[i]);
                    }
                }
            }
            return result;
        }


        public Vector S()
        {
            return mS.Clone();
        }

        public Matrix Solve(Matrix input)
        {
            if (!mComputeVectors)
            {
                throw new InvalidMatrixOperationException(Resources.SingularVectorsNotComputed);
            }
            if (!Converged())
            {
                throw new ConvergenceFailedException();
            }
            Matrix result = mU.CreateMatrix(mColumns, input.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Matrix input, Matrix result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (input.Columns != result.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (!mComputeVectors)
            {
                throw new InvalidMatrixOperationException(Resources.SingularVectorsNotComputed);
            }

            if (!Converged())
            {
                throw new ConvergenceFailedException();
            }

            if (mU.Rows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mV.Columns != result.Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            int mn = System.Math.Min(mRows, mColumns);
            int bn = input.Columns;

            double[] tmp = new double[mColumns];

            for (int k = 0; k < bn; k++)
            {
                for (int j = 0; j < mColumns; j++)
                {
                    double value = 0;
                    if (j < mn)
                    {
                        for (int i = 0; i < mRows; i++)
                        {
                            value += mU.ValueAt(i, j)*input.ValueAt(i, k);
                        }
                        value /= mS[j];
                    }
                    tmp[j] = value;
                }

                for (int j = 0; j < mColumns; j++)
                {
                    double value = 0;
                    for (int i = 0; i < mColumns; i++)
                    {
                        value += mV.ValueAt(i, j)*tmp[i];
                    }
                    result[j, k] = value;
                }
            }
        }

        public Vector Solve(Vector input)
        {
            if (!mComputeVectors)
            {
                throw new InvalidMatrixOperationException(Resources.SingularVectorsNotComputed);
            }
            if (!Converged())
            {
                throw new ConvergenceFailedException();
            }
            Vector result = mU.CreateVector(mColumns);
            Solve(input, result);
            return result;
        }

        public void Solve(Vector input, Vector result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (!mComputeVectors)
            {
                throw new InvalidMatrixOperationException(Resources.SingularVectorsNotComputed);
            }

            if (!Converged())
            {
                throw new ConvergenceFailedException();
            }

            if (mU.Rows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mV.Columns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            int mn = System.Math.Min(mRows, mColumns);
            double[] tmp = new double[mColumns];
            double value;
            for (int j = 0; j < mColumns; j++)
            {
                value = 0;
                if (j < mn)
                {
                    for (int i = 0; i < mRows; i++)
                    {
                        value += mU.ValueAt(i, j)*input[i];
                    }
                    value /= mS[j];
                }
                tmp[j] = value;
            }

            for (int j = 0; j < mColumns; j++)
            {
                value = 0;
                for (int i = 0; i < mColumns; i++)
                {
                    value += mV.ValueAt(i, j)*tmp[i];
                }
                result[j] = value;
            }
        }

        private void Compute()
        {
            if (!mComputed)
            {
                DoCompute();
                mComputed = true;
            }
        }

        protected abstract void DoCompute();
    }
}